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In quantum simulation, many-body phenomena are probed in controllable quantum systems. Re¬ 
cently, simulation of Bose-Hubbard Hamiltonians using cold atoms revealed previously hidden local 
correlations. However, fermionic many-body Hubbard phenomena such as unconventional super¬ 
conductivity and spin liquids are more difficult to simulate using cold atoms. To date the required 
single-site measurements and cooling remain problematic, while only ensemble measurements have 
been achieved. Here we simulate a two-site Hubbard Hamiltonian at low effective temperatures 
with single-site resolution using subsurface dopants in silicon. We measure quasiparticle tunneling 
maps of spin-resolved states with atomic resolution, finding interference processes from which the 
entanglement entropy and Hubbard interactions are quantified. Entanglement, determined by spin 
and orbital degrees of freedom, increases with increasing covalent bond length. We find separation- 
tunable Hubbard interaction strengths that are suitable for simulating strongly correlated phenom¬ 
ena in larger arrays of dopants, establishing dopants as a platform for quantum simulation of the 
Hubbard model. 


INTRODUCTION 

Quantum simulation offers a means to probe many- 
body physics that cannot be simulated efficiently by clas¬ 
sical computers, using controllable quantum systems to 
physically realize a desired many-body Hamiltonian 1 3 . 
In the analog approach to quantum simulation exem¬ 
plified by cold atoms in optical lattices 4 the simu¬ 
lator’s Hamiltonian maps to the desired Hamiltonian. 
Compared to digital quantum simulation, realized via 
complex sequences of gate operationeP® analog quan¬ 
tum simulation is usually carried out with simpler build¬ 
ing blocks. For example, the Heisenberg and Hubbard 
Hamiltonians of great interest in many-body physics are 
directly synthesized by cold atoms in optical lattices 23 . 
Although of immense interest and proposed long ag(P, 
analog simulation of fermionic Hubbard systems has 
proven to be very challenging 2 3 . The anticipated regime 
of the intensely debated spin liquid, unconventional su¬ 
perconductivity, and pseudogapPC 11 has yet to be ac¬ 
cessed even for cold atoms. Here, the required low 
temperature T < t/30 is probl emat ic due to the weak 
tunnel-coupling t of cold atomd^ 12 . Moreover, experi¬ 
mentally resolving individual lattice sites, crucial else¬ 
where in Bose-Hubbard simulation®, remains very chal¬ 
lenging in quantum simulation of the Hubbard model. 

Here, we perform atomic resolution measurements re¬ 
solving spin-spin interactions of individual dopants, re¬ 
alizing an analog quantum simulation of a two-site Hub¬ 
bard system. We demonstrate the much desired combi¬ 
nation of low effective temperatures, single-site spatial 
resolution, and non-perturbative interaction strengths of 
great importance in condensed matter 9 11 . The dopants’ 
physical Hamiltonian PL S i m , determined at the time of 
fabricatiorP, maps to an effective Hubbard Hamiltonian 


Usys = Ei^>(^ C L C J> + h - C 0 + 'Ei,a Un it n iP where 
U is the on-site Coulomb repulsion, cj a ( Ci a ) creates (de¬ 
stroys) a fermion at lattice site i with spin cr, = ct^c^ 
is the number operator, and h.c. is the Hermitian con¬ 
jugate. Here, it is desirable to achieve non-perturbative 
(intermediate) interaction strengths U/t associated with 
quantum fluctuations and emergent phenomena 9 ®®, z.e., 
beyond perturbative Heisenberg interactions (large U/t) 
realized in photon-based 13 and ion-basec£ 14 l simulations, 
and magnetic ions on metal surface! 15 . We focus on 
the system ground state, prepared by relaxation upon 
cooling®, rather than system dynamics. 

Because the states of our artificial Hubbard system 
are coupled and interacting, tunneling spectroscopy lo¬ 
cally probes the spectral function. The spectral function 
is of key interest in many-body physics because it pro¬ 
vides rich information on interactions 16 17i , and is highly 
sought after in future “cold-atom tunneling microscope” 
experiments^! For our few-body system, the local spec¬ 
tral function describes the quasi-particle wavefunction 
(QPWF)P^ 22 I and the discrete coupled-spin spectrum of 
the dopants. We find that interference of atomic orbitals 
directly contained in the QPWF allows us to quantify 
the electron-electron correlations and the entanglement 
entropy. The entanglement entropy is a fundamental con¬ 
cept for correlated many-body phases 23 ® that has thus 
far evaded measurement for fermions. In the counterin¬ 
tuitive regime of our experiments, entanglement entropy 
increases as the valence bond is stretched, as Coulomb 
interactions overcome quantum tunneling. In our sys¬ 
tem, the entanglement entropy is directly related to the 
Hubbard interactions U/t , and we find that U/t is tun¬ 
able with dopant separation, increasing from 4 —» 14 
for d/ctB = 2.2 —)> 3.7, where a b = 1.3 nm is the ef¬ 
fective Bohr radius. This range, of interest to simulate 
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FIG. 1. Spatially resolving coupled-spin states A. 

Atomic resolution single hole tunneling probes the interact¬ 
ing states of two coupled acceptor dopants (r ou t= tunnel 
rate to tip, r ou t C Fn = tunnel rate from reservoir) The 
inter-acceptor coupling t obeys t hT m . dl/dU measures 
the interacting states’ quasi-particle wavefunction (QPWF), 
contains interference processes from which we obtain two- 
body wavefunction amplitudes, and determine the entangle¬ 
ment entropy and effective Hubbard interactions. B. Accep¬ 
tor pair (double-protrusion) in topography at U — +1.8 V 
and I = 300 pA (top), and spectrally and spatially resolved 
dl/dU taken at a bias U = +2.0 V where topography is flat 
apart from atomic corrugation (bottom). Valence band (VB), 
2-hole ground state and 2-hole excited states are indicated. C. 
Effective energy diagram of sequential hole tunneling through 
2-hole ground and excited state of coupled acceptors. 


unconventional superconductivity and spin liquidd^J] i s 
realized here due to the large Bohr radii of the hydro- 
genic states. The semiconductor host allow s for electro¬ 
static control of the chemical potential 2 ^!!, desirable to 
dynamically control filling-factor 9 11 but not possible for 
ions on metal surfaced^. 


RESULTS 


Spectroscopy of coupled-spin system Subsurface 
boron acceptors i n silic on were identified at 4.2 K as indi¬ 
vidual protrusions 29 30 (density ^ 10 11 cm -2 ) in constant 
current images due to resonant tunneling at a sample bias 
U = +1.6 V, and due to the acceptor ion’s influence on 
the valence density of states at U = —1.5 V. The sample 
was prepared by ultra high vacuum flash annealing at 
1200 °C and hydrogen termination. The observed sub¬ 
surface acceptors had typical depthJ^SEo] < 3 nm? an d 
correspondingly, a volume density > 25 times less than 
the bulk doping, 8 x 10 18 cm -3 . Pairs of nearby acceptors 
with d < 5 nm were also found, with a smaller density 
~ 10 9 cm -2 . 

The spectrum and spatial tunneling probability of the 
coupled acceptors were investigated at T = 4.2 K via 
single-hole tunneling from a res ervoir in the substrate to 
the dopant pair, to the tipl 29 | 3Q | (Fig. [+)• For the dopant 
pair in Fig. (top), dl/dU measured along the inter¬ 
dopant axis (Fig. [lj3, bottom) contains a peaks for each 
state entering the bias window, at U ^ 0.2 V, 0.45 V, 0.55 
V and 0.8 V. Consistent with our single-acceptor 29 and 
single-donor 31 measurements near flat-band bias condi¬ 
tions, the bias for each peak in the spectrum (Fig. [l]3, 
bottom) is independent of tip position. This rules out 
distortion of our quantum state images by inhomogenous 
tip-induced potential ^ 2 observed in other multi-dopant 
systems 33 . These results can be attributed to weak elec¬ 
trostatic control by the tip (Fig. FT) and the states’ prox¬ 
imity to flat-band 2 ^®], though a large tip radius may also 
play a role. 

The spectral and spatially resolved measurements 
(Fig. [l]3) directly demonstrate that the holes are 
interacting, as follows. First, two peaks centred on 
dopant ions A or B are resolved in real space (Fig. [+)• 
Second, energy differences between the peaks resolved 
in real space are smaller than the ^ 350 /ieV thermal 
resolution. However, for orbitals at the same energy 
to not interact, their overlap must vanish. Since the 
measured orbitals have a strong overlap, the sites are 
tunnel coupled, irrespective of the details of the tun¬ 
neling current profile. The number of states observed, 
their energy differences, and their energies relative to 
the Fermi energy confirm that they observed states are 
two-hole states (Fig. [ 6 ] and [7]). 

Correlations, Entanglement, and Hubbard Inter¬ 
actions The ground state of a Hubbard model with non- 
perturbative interactions is governed by H in Fig. [ 2 ]A in 
the subspace of |t;|) = c^c^| 0 ), ||;t) = I 0 )’ 

I'M-;) = ++ | 0 ) and |; U) = +<+ | 0 ), where + ere- 
ates a localized electron on site i G {A, B} with spin a G 
{t, and |0) is the vacuum state. The ground state is a 
superposition |^ s ) = 7c(|tU) - 14-51)) +7i(|U;} + |; tl», 
where + c (+i) is the probability amplitude for a cova¬ 
lent (ionic) configuration (Fig. [2^). Rewriting the state 
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FIG. 2. Hubbard Interactions and Entanglement En¬ 
tropy A. Two-site Hubbard Hamiltonian in the subspace of 
the ground state, with tunnel coupling t hybridizing singly 
and doubly occupied configurations, for sites A (red orbital) 
and B (blue orbital). B. Dependence of probability ampli¬ 
tudes on interactions hi/t: y c (green dashed) and 71 (green 
solid) for configurations (|t;>l) — \Ut)) and (\tU) + |;t^))> 
and 7ee and 7 00 for |e^e^) and |o^o^} respectively. C. Entan¬ 
glement entropy S increases with increasing Hubbard interac¬ 
tions U/t. This occurs because of localization of red and blue 
orbitals associated with spins in the singlet, as illustrated in 
the insets. 


in a basis of even and odd orbitals, |\hs) = 7 ee l e t e t) — 
7 oo |o^o^), where 7 ee ( 700 ) is the probability amplitude 
of the “even/even” (“odd/odd”) configuration. 

In limit of small tunnel couplings (U/t -A 00, Fig. | 2 ^) 
the Hubbard system may be described by perturba¬ 
tive Heisenberg spin interactions. For vanishing U/t, 
the ground state is a Heitler-London singlet of local- 
ized spins, |^ s ) = 2 _ 1 / 2 (|T; 4 ) - 14-5T», with no contri¬ 
butions from |H;) and Due to vanishing wave- 

function overlap the electrons can be associa ted with 
sites A and B (they are distinguishable ^ 34 * 35 * ) , and the 
spin at site A depends on the spin at site B as for a 
maximally entangled Bell state. In the limit of vanish¬ 
ing interactions ( U/t -T 0, Fig. [2^3) corresponding to 
a tight-binding approximation, the spins delocalize and 

l*s> = Kit; 4-) - I4-;t» + §(lt4-;> + l;t4-». in a molec- 

ular orbital (MO) basis, the ground state |Ts) = |e^), 
which is a single Slater determinant. Although this state 
is a singlet (one spin up, one spin down) due to funda¬ 
mental indistinguishability, the electrons can be ascribed 
independent properties because th ey occu py the same or¬ 
bital, and the state is nncorrelate d 23 * 34 *^. 

In the regime of intermediate U/t where tunneling and 
Coulomb interactions compete non-pertnrbativel ^ 2 * 3 * 9 * 11 *, 
tunnelling hybridizes the doubly-occupied configurations 
|ti;) and |; tl) into the ground state, such that the 
particles lose their individual identities. Here, the 
von Neumann entanglement entropy quantifies gen¬ 


uine entanglement (inter-dependency of properties), 
distinguishing it from exchange-correlations due to 
indistinguishabilit}^ 3 ! 26 ! 35 ] Employing the convention— 
S = 0 ( 1 ) for zero (maximal) entanglement, S = 
-| 7 ee| 2 log 2 | 7 ee | 2 - | 7 oo| 2 log 2 | 7 oo | 2 increases as U/t in¬ 
creases and coherent localization occurs (Fig. [2p) , satu- 
rating at value of 1 . 

We now discuss the spatial tunneling maps of the two- 
hole ground states for different inter-acceptor distances. 
Obtained by integrating the lowest voltage d I /d U peak, 
the maps are shown in Fig. [3^, [3^3, and HP for dis¬ 
tances d/a& = 2.2,2.7 and 3.5 (ap = 1.3 nm) hav¬ 
ing orientations ± 2 ° from ( 110 ), 8 ± 2 ° from ( 100 ) and 
3 zb 2° from (110), respectively. The multi-nm spatial 
extent of the states reflects the extended wave-like na¬ 
ture of the acceptor-bound holes, owing to their shal¬ 
low en ergy levels, which contrasts Mn ions on GaAs 
surface!®*, magnetic ions on metald®*, and Si( 001 ):H 
dangling bonds 38 . Consequently, their envelopes are 
amenable to ef fective-ma ss analysis with lattice frequen¬ 
cies filtered out 19 2 ®SI 39 [ Consistent with measurements 
of singl e acc eptors at similar depths on resonance at 
flatbancl 2 ^^, the states have predominantly s-like en¬ 
velopes with slight extension along [ 110 ] directions, as 
expected when symmetry is not strongly perturbed by 
the surface. Depths of the d/ap = 2.7 and d/ap = 3.5 
pairs were estimated to be ~ 0.9 nm, and for d/ap =2.2, 
~ 0.6 nm (see Fig. [ 8 |. 

We employed full-configuration interaction calcula¬ 
tions of the singlet ground-state |\hs) to confirm that 
Coulomb correlations of coupled acceptors influence the 
ground state in a way that mimics the 5=1/2 Hubbard 
model. In particular, for d/ap ^ 2, |T$) is predomi¬ 
nantly composed of c £ 3 / 2 c e -3/2 1 ^)’ a s i n gl e ^ of two even 
±“3/2” spin MOs. With increasing d, interactions en¬ 
hance the probability amplitude of the c/ 3 / 2 c o —3/2 10 ) 
singlet with two odd orbitals, analogous to the Hub¬ 
bard Hamiltonian (Fig. [2b). The spins ±“3/2” are pre¬ 
dominantly composed of|3/2, ±3/2) valence band Bloch 
states. In particular, the low-lying ±“1/2” spin excita¬ 
tions of each acceptor®!, w hich are predominantly com¬ 
posed of |3/2,±l/2) Bloch states, do not qualitatively 
change the description. We also note that for d/ap > 2 , 
the MOs are essentially linear combinations atomic or¬ 
bitals having the effective Bohr radii of single acceptors. 

Single-hole tunneling transport through our coupled 
dopant system locally probes the spectral quasi-particle 
wavefunctior^^. When r out r in (Fig. [lpA) , the 
single-hole tunneling rate is essentially governed by r out , 
the tunnel-out rate 31 . In the present case, single-hole 
tunneling from the two-hole system to a single-hole fi¬ 
nal state I/) = cjr|0) (Fig. m) contributes rl ut (r) = 
|(/|^(r)|* s )| 2 , where (/|T(r)|* s ) is the QPWF, T(r) = 
(f)j(r)cj is the field operator, c| creates a single-hole 
MO eigenstate <fij (r) of the system^, and the total tunnel 
rate is T(r) -£/ r o Ut (r)- 

From the QPWF description of coupled dopants, we 
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FIG. 3. Resolving Interference Processes in Quasi-Particle Wave Function A. Experimentally measured, normalized 
tunneling probability T oc dl/dU to tip, for d = 2 . 2 as ground state. Arrows denote 110 crystal directions. B. Same as (A), for 
d — 2.7ae. C. Same as (A), for d = 3.5as. D. Normalized experimental line profile (coloured squares) of T(x) for d = 2.2ub and 
least-squares fit (coloured line) to QPWF correlated singlet model. Lower and upper grey lines are line profiles of maximally 
and minimally correlated states, obtained from least square fits. The maximally correlated state deviates from the mean of 
| (J) e (r) | 2 and |0 o (r)| 2 because of the different normalization coefficients of even and odd linear combinations. E. Same as (D), 
for d m 2.7as- F. Same as (D), for d = 3.5as- Scale: 1 nm. 


obtain a spatial tunneling probability T(r, | 7 ee |, | 7 oo|) oc 
|7ee| 2 |0e(r )| 2 + | 7 oo| 2 |4>(r )| 2 for the ground state. Here, 

17ee1 2 an d (|7 00 | 2 ) contain constructive (destructive) 
interference corresponding to even (odd) linear com¬ 
binations of atomic orbitals </> e ( r i) ( 0 o ( r i)) (note: 

17ee| 2 + | 7 oo | 2 = 1)- To obtain | 7 00 | 2 , data were 

fit to T(r, 17 ee 15 | 7 oo|)? assuming linear combinations of 
parametrized s-like atomic orbitals for </> e (r) and 0 o (r) 
appropriate for subsurface acceptors. The QPWF and 
atomic orbitals are described in Figs. [9| [To| and [TTj 

The least-squares fits in Figs. HE Fp, [3p (colored 
lines) of T(r, 1 7 ee |, 17 00 1) are in good agreement with data 
(squares), for d/a-g = 2.2, 2.7 and 3.5. For comparison 
with the data, grey curves are shown for both the un¬ 
correlated (maximally correlated) state with | 7 00 | = 0 
(|7oo|/|7ee| = 1) in Fig. |^D, [3p and |T We note that 
all three separations exhibit interaction effects at the 
midpoint of the ions, where the quantum interference is 
strongest. We obtain 17 00 1 2 = 0.12 db 0.06, 0.23 4= 0.07, 
and 0.39T0.08 for d/ap = 2 . 2 , 2.7 and 3.5. Data taken at 
higher tip heights gave identical results to within experi¬ 
mental errors (see Fig. [I 2 | and 13), independently verify¬ 
ing that the tip does not influence our results. 


The Coulomb correlations, embodied both in C = 
2|7oo| 2 (Fig- [fj^) and the entanglement entropy S = 
-|7ee| 2 log 2 l7ee| 2 - l7oo| 2 fog 2 l7oo| 2 (Fig. [4p) , could be 
evaluated directly from the fit, and both increase with 
increasing d. The one-to-one mapping from S to U/t 
(Fig. IT) was used to determine the effective Hubbard 
interactions from the entanglement entropy in Fig. [4|3. 
We obtain U/t « 3.5, 6.4, and 14, for d/ap = 2 . 2 , 2.7, 



Separation (a g ) 


FIG. 4. Entanglement Entropy and Hubbard inter¬ 
actions A. Quantum correlations C vs. d. Theory predic¬ 
tions are shown for coupled acceptors with ( 110 ) orientations 
(red line) and ( 100 ) (blue line), alongside scaled H 2 (dashed 
black line). Predicted localization is suppressed (enhanced) 
along ( 110 ) (( 100 )) relative to molecular hydrogen (H 2 ), due 
to valence band anisotropy, which enhances (suppresses) t. B. 
Same as (A), for the entanglement entropy S. C. Experimen¬ 
tally estimated Hubbard interactions. Error bars denote 95 
% confidence intervals. 
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and 3.5 respectively (Fig.[4p), which increase as the tun¬ 
nel coupling decreases. 

We conclude the analysis of the QPWFs with some 
critical remarks on correlations extracted from our 
fitting model, recalling that the large spatial overlap 
of the spectrally overlapping acceptor-bound holes 
directly shows their states are tunnel coupled. First, 
the Coulomb correlations have a systematic effect on 
interference in the QPWF such that the least-squares 
error is significantly worse if | 7 00 | 2 is forced to zero in 
the fitting model (Table |T|. Second, if applied to very 
far apart dopants where the ground state can still be 
resolved, our fitting model would not give a spurious 
result that the two dopants are highly correlated. 
This follows because the difference between |</> e (r)| 2 
and |0 o (r)| 2 , which reflects the interference of atomic 
orbitals and is used to detect correlations, tends to 
zero as d/aB increases. Data (Fig. [± -C) presented 
here are for coupled dopants that we found to be (i) 
well isolated from other dopants or dangling bonds, 
and (ii) at identical depths, as evidenced by the spatial 
extent and brightness of the atomic orbitals. When the 
latter is not satisfied, the atomic levels can be detuned, 
introducing more parameters to the fit. 

Comparison with theory These experimental re¬ 
sults obey the trends predicted by our theory calculations 
for the spin-orbit coupled valence band. Predictions in 
Fig. [4pV and [4b for displacements along (100) (blue solid 
line) and ( 110 ) (red solid line) both show increasing cor¬ 
relations and entanglement with increasing dopant sep¬ 
aration. Moreover, we find that the observed and pre¬ 
dicted entanglement entropy qualitatively reproduce a 
single-band model (Fig. [4|4, [4^3, dashed lines). This re¬ 
sult implies that inter-hole Hubbard interactions follow 
an essentially hydrogenic trend with atomic separation, 
even for non-perturbative interactions U/t — 4 ^ 14. 

The hydrogenic nature of S and U/t persists in spite 
of the ±“ 1 / 2 ” spin excited states of a single acceptors. 
Such ±“ 1 / 2 ” single-acceptor excited states states are 
found nominally A ^ 1 — 2 meV above the ±“3/2” 
spin ground state due to inversion symmetry breaking 
at the interfacd^. Although t > A, S and U/t remain 
hydrogenic in our calculations because the “ 1 / 2 ” spin 
excited state has an s-like envelope whose spatial extent 
is similar to (1) the s-like ±“3/2” ground state and (2) 
the scaled hydrogenic ground state. Otherwise, single 
particle ±“ 1 / 2 ” states would hybridize stronger than 
single particle ±“3/2” states, form the 2-hole singlet at 
smaller separations, and localize more slowly relative 
to molecular hydrogen with increasing d. Furthermore, 
the polarization of the ±“3/2” and ±“1/2” states into 
|3/2, ±3/2) and |3/2, ±1/2) components, respectively, 
limits the mixing of ±“ 1 / 2 ” states into the ground state. 

Spin excited states and effective temperature 

Finally, we discuss the observed excited states, which 
confirm that the inter-acceptor tunnel-coupling domi- 
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FIG. 5. Coupled-spin excitation spectrum A. Measured 
energy of first excited state relative to ground state. B. 
Schematic level diagram of coupled acceptors, reflecting the¬ 
ory calculations, as a function of inter-acceptor distance d/a-Q. 
Singlets |£ mj ) and triplets |T mj ) are present for interactions 
between two holes of mj = ±“3/2” spin (orange) and two 
holes of mnj — ±“ 1 / 2 ” (black) spin. States | Q3/2P/2) and 
\Q 3/2,1/2) are se ^ s °f f° ur closely spaced levels (green) with 
one “3/2“ spin hole, and one “1/2” spin hole. Error bars 
denote 95 % confidence intervals. 


nates thermal and tunnel-coupling effects of the reservoir. 
The energies of the states were determined by fitting the 
single-hole transport lineshapes 40 of the coupled accep¬ 
tors (Fig. [6] and [7|). For the first excited state we found 
5.2 ± 0.6 meV and 1.2 ± 0.2 meV for d/a b = 2.2 and 
3.5 respectively (^ (110) orientation), and 1.6±0.7 meV 
for d/a-Q = 2.7 (~ (100) orientation). Shown in Fig. [ 5 JA, 
these energies are too small to add another hole, which 
would require « 50 meV for an acceptor in bulk silicon. 
However, the energies agree well with our predictions for 
two-hole excited states of coupled hole spins ±“3/2” and 
±“1/2”, i.e., 8.5 meV and 1.5 meV for d = 2.2ae and 
d = 3.5<2b ((HO) orientation), and 2.0 meV ((100) orien¬ 
tation). Here we note that some of the predicted coupled- 
spin excited states (Fig. IT) are unconventional: a sin¬ 
glet | Sjjij ) and triplet |T mj ) of two “3/2” holes (orange 
lines) and two “1/2” holes (black lines), where l-S' 3 / 2 ) is 
the ground state for all separations. More subtly, two 
manifolds \Qy 2 1/2 )> 1 * 33 / 2 , 1 / 2)5 i — 1 - - . 4, containing 
four states are predicted (green lines), where one ±“3/2” 
spin level and one ±“1/2” spin level is occupied. For 
d/aB = 2.2 and 2.7 (d/a-Q = 3.5), the measured energies 
are in better agreement with predictions for \Q\/ 2 1 / 2 ) 
(|T 3 / 2 » excitations. 

The inter-acceptor tunnel couplings t (ratios t/T) were 
estimated to be 12 meV (30), 7 meV (20), and 3.5 meV 
(10) for d/aB = 2.2, 2.7 and 3.5 respectively, at T = 4.2 
K. Such couplings t exceed the reservoir coupling Ti n (Ta¬ 
ble 0 to the substrate by more than 50X. Combined 
with bias U ~ 0.2 — 0.3 V needed to bring the level 
into resonance, this rules out coherent interactions with 
substrate and tip reservoird^. Note that the measured 
energy splittings imply small thermal excited-state popu- 
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DISCUSSION 


AUTHOR CONTRIBUTIONS 


We performed atomic resolution measurements resolv¬ 
ing spin-spin interactions of interacting dopants, real¬ 
izing quantum simulation of a two-site Hubbard sys¬ 
tem. Analyzing these local measurements of the spectral 
function^, we find increasing Coulomb correlations and 
entanglement entropy as the system is “stretched” 1 - 3 33 42 ^ 
in the regime of non-perturbative interaction strengths 
U/t. Our experiment is the first to combine low effective 
temperatures t/T ~ 30 at 4.2 K and s ingle- site mea¬ 
surement resolution, considered essential 3 5 12 to simu¬ 
late emergent Hubbard phenomena 911 . Lower effective 
temperatures t/T ~ 420 are possible at T = 0.3 K. For 
example, 4x4 Hubbard lattices with U/t = 4 7 and 

t/T ~ 40 have recently been associated with both the 
pairing state and pseudogap in systems exhibiting un¬ 
conventional superconductivity] 

The approach generalizes to donors, which can be 
placed in silicon with atomic-scale precision 27 and spa¬ 
tially measured in-situ after epitaxial encapsulation 43 44 . 
In contrast to disordered systemd^, atomically engi¬ 
neered dopant lattices will require weak coupling to 
a reservoir, displaced either vertically as demonstrated 
herein, or a lateralljJ^. Strain could be used to fur¬ 
ther enhance the splitting between light and heavy h oles, 
or suppress valley interference processes of electrong i 31 * 46 [ 
Interestingly, open Hubb ard sy stems which may exhibit 
unusual Kondo behaviou j 47 1 48 * could also be studied by 
this method. The demonstrated measurement of spectral 
functions could be used to directly determine excitation 
spectra, evaluate correlation function^, or obtain quasi¬ 
particle interference spectra 1 -^, all of which contain rich 
information about many-body states, including charge¬ 
ordering effects. We envision in-situ control of filling 
factor^], using a back-gate or patterned side-gate 27 . 
These capabilities will allow for quantum simulation of 
chains, ladders, or latticed 9 at low effective temper¬ 
atures, having interactions that are engineered atom-by¬ 
atom. 


Experiments were conceived by J.S, J.A.M, and S.R. 
J.S. carried out the experiments and analysis, with input 
from J.A.M., R.R., L.C.L.H, and S.R. Theory modeling 
was carried out by J.S., J.A.M., R.R., and L.C.L.H and 
S.R., with input from all authors. J.S. and S.R. wrote 
the manuscript with input from all authors. 
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METHODS 


Sample Preparation Samples were prepared by flash 
annealing a boron doped (p « 10 19 cm -3 ) silicon wafer 
at ~ 1200 °C in UHV followed by slow cooling at a 
rate 1 °C-min _1 to 340 °C. Then, hydrogen passivation 
was carried out ^ 340 °C for ten minutes by thermally 
cracking H 2 gas at a pressure Ph 2 = 5 x 10“ 7 mbar. 


Measurements Atomic resolution single-hole tunneling 
spectroscopy was performed at 4.2 K using an ultra-high 
vacuum Omicron low temperature scanning tunneling 
microscope (LT-STM). Current / was measured as a 
function of sample bias U and dl/dU was obtained by 
numerical differentiation. Details for the analysis of the 
data are provided in Fig s. [6} [7| [8j [lOj |ll|12[ and [~L3| and 
Appendices | A| B| D| and \M 


Theory Theory calculations of interacting states were 
carried out using the configuration interaction approach, 
in the Luttinger-Kohn representation including a realis¬ 
tic description of the heavy-hole (J = 3/2,|raj| = 3/2), 
light-hole (J = 3/2,|raj| = 1/2), and split-off hole 
(J = 1/2,|raj| = 1/2) degrees of freedom. Details for 
the theory are provided in Fig. [9] and Appendices | C| F 
andl~Gl 


ACKNOWLEDGEMENTS 

We thank H. Wiseman, M. A. Eriksson, M. S. Fuhrer, 
O. Sushkov, D. Culcer, J.-S. Caux, B. Reulet, G. 
Sawatzky, J. Folk, F. Remade, M. Klymenko and B. 
Voisin for helpful discussions. This work was supported 
by the European Commission Future and Emerging Tech¬ 
nologies Proactive Project MULTI (317707), the ARC 
Centre of Excellence for Quantum Computation and 
Communication Technology (CE110001027), and in part 
by the U.S. Army Research Office (W911NF-08-1-0527) 
and ARC Discovery Project (DP120101825). S.R. ac- 


SUPPLEMENTARY INFORMATION 
A. Tunneling spectroscopy 

In this section we discuss the extraction of the energies 
of the states observed in dl/dU tunneling spectra. We 
find that the the number of observed states, the state’s 
relative energies to eachother, and their energies relative 
to the sample’s Fermi energy, are only compatible with 
two interacting holes on pairs of acceptors. We further 
estimate the total tunnel coupling T^ to the tip and sam¬ 
ple reservoirs for sequential hole transport, based on the 
dl/dU lineshape, and show that hT^ is much less than 
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FIG. 6. A. dl/dU vs. U for the d/dB — 2.2, above (black 
squares) and away from (orange squares) the coupled accep¬ 
tors. Least square fits (green line) considering sequential hole 
transport through several peaks (green dashed lines). B. Same 
as A on logarithmic scale. C. Tunneling of electrons from va¬ 
lence band (v.b.) to the tip. D. Tunneling of electrons from 
the tip to the conduction band (c.b.). Schematic of hole tun¬ 
neling at positive sample bias U, from valence impurity band 
reservoir, to acceptor pair, to tip. Ground and excited states 
for two-hole occupation are shown. Electrons (holes) in the 
hole reservoir and STM tip are shaded blue (white). 


the inter-acceptor tunnel coupling t for separations d con¬ 
sidered in the main text. 

For the acceptor pair with d/a-g = 2.2 (Fig. [3K, main 
text), measured dl/dU is plotted in linear scale (Fig. [6|\) 
and log scale (Fig. §3). Away from the acceptor pair 
(Fig.[6]A/B, orange squares), direct tunneling of electrons 
from the valence band to the tip (holes from the tip to 
the valence band) is obtained for U < 0 V as illustrated 
in Fig. [6p, while direct tunneling of electrons from the 
tip to the conduction band occurs for U > 1.15 V as 
illustrated in Fig. J6p. Above the pair (Fig. [6]A/B, black 
squares), we observe a well-isolated peak at U ~ 0.3 V, 
a collection of two strong (closely spaced) peaks at 0.6 V 
and 0.75 V, and fourth peak at U ^ 1.0 V. 

Plotted on a logarithmic scale, both the current 
(not shown) and numerically differentiated dl/dU data 
(Fig. [6^), black squares) grow exponentially with in¬ 
creasing V for the first peak before reaching a local 
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FIG. 7. Extracted energies Ei — E\ of observed states relative 
to the lowest energy state (energy Ffi), for d/dB — 2.2, d/dB — 
2.7, and d/dB — 3.5. For comparison, the energies of the 
predominantly ±“3/2” and ±“1/2” spin Kramers doublets 
of a lone acceptor, ~ 1 nm from the Si[001]:H interface, are 
shown on the right hand side. 


maximum. As suggested by this observation 29 ^51, we fit 
the first peak of the differential conductance to a state 
probed by sequential tunneling with broadening deter¬ 
mined by the T = 4.2 K temperature of the reservoir, 
cosh ~ 2 (ea(U — U{)/2kBT), where e = 1.602 x 10 -19 C 
is the elementary charge, = 1.381 x 10 -23 J/K is 
the Boltzmann constant, U\ is the peak voltage, and a 
is the well-known lever arm from sequential transport— 
describing the variation in the energy level A E #= eaAU 
due to a change in applied bias Af7, as schematically 
illustrated in Fig. [6]E. 

The least-squares fit is in good agreement with the 
measured d I / dU for both the full-width at half maximum 
A U = 3 .dhsT/ea and the exponential decay of the tail 
exp (ea(U — Uij/ksT). A peak voltage U\ = 334.0 ± 0.9 
meV and a lever arm a = 0.0144 ± 0.0006 were ob¬ 
tained for tunneling of holes from the reservoir, to the 
acceptor pair, to the tip. The value for a is similar 
to those found for single-hole transport through individ¬ 
ual acceptors in flashed p-type silicorl 29 and single ar¬ 
senic donors in flashed n-type silicon 31 . This excellent 
agreement shows that broadening of the single-hole trans¬ 
port peaks is due to the temperature associated with the 
Fermi-Dirac distribution in the hole reservoir. 

The small value of a means that the “broadened” line- 
shape in Fig. Hk and Fig. [6j3 is due to weak electro¬ 
static control of the sample potential by the STM tip, 
as expected due to strong screening by the carrier reser¬ 
voir. This is further confirmed by the negligible “spectral 
shift’®, that is, the similarity of the apparent gap « 1.15 
eV in dl/dU measurements of Fig. [6jA/B to the actual 
1.14 eV band gap for heavily-doped p-type silicon at low 
temperature 52 . 




































1. Energy spectrum 


We extracted the lever arm a and energies of states 
Ei~E\ = ea(Ui — Lh), by least-squares fit ting of dl/dU 
to a sum of single-hole tunneling peaks^IHSi [40 53 ] Fits 
to dl/dU = Ai(z)(cosh~ 2 (ea(U - Ui)/2k B T )) are 
found in Fig. pA and pfe (solid green lines) along with 
individual peaks (dashed green lines), on linear and log 
scales, respectively. Extracted energies Ei~E\ — ea(Ui — 
Ui) of the excited states are plotted in Fig. [7| relative to 
the energy Ei of the two-hole ground state, for d/a b = 
2.2, d/aB = 2.7 and d/aB = 3.5. 

For all d the extracted energy splittings between the 
states (Fig. [7]) are too small to be associated with single¬ 
hole occupied states (A^T states) of coupled acceptors. 
For d/aB = 2.2, two energetically nearby states were 
found at 5 and 6.5 meV, while a higher energy state was 
found at ~ 3.5 meV above those. Here, the observed 
splittings are incompatible with the “Is” hybridization 
energy 2 1 ~ 25 meV between tunnel-hybridized single¬ 
hole states estimated by numerics (Section [C]). Similar 
arguments hold for the d/aB = 2.7 (d/aB = 3.5), where 
the excited state splittings « 2 meV (« 1 meV) are in¬ 
compatible with the single-hole coupling energy scale of 
2 1 ~ 14 meV ( 2 1 ~ 7 meV). 

Even more importantly, it follows from the electro¬ 
static arguments given below that the coupled-acceptor 
peaks observed in experiments are too close to the hole 
reservoir’s chemical potential to correspond to ionizations 
of the very “deep” single-hole (A^ -1 ) states into zero-hole 
(A^ 2 ) states. Rather, they are only compatible with ion¬ 
izations of the two-hole (A 9 ) ground and excited states 
of acceptor pairs into one-hole (A^ 1 ) states, which are 
energetically very similar to the ionization energy of a 
neutral (A 0 ) acceptor. For d/aB = 2.2 —» 3.5 we found 
ionization transitions A ^ 1 A 2 2 for the last hole (Sec¬ 
tion [C| require more than 100 meV energy. This exceeds 
the Fermi energy Ep — Ey 50 meV of the reservoir ( Ey 
is the valence band edge) by more than 50 meV. Now, 
for U = 0 V, the tip’s contact potential bends states by 
an amount — eaUpB = — ea(W — <F t i p ) relative to the 
sample reservoir, where W = 5.2 eV is the work function 
of degenerately doped p-type reservoir, and 4> tip is tung¬ 
sten tip workfunction. Let us assume a work function 
3>tip = 4.8 ±0.3 eV typical of a tungsten tipPP. Then, the 
states bend down by 5.8 d= 4.2 meV at zero bias, much 
smaller than the value > 50 meV required to depopulate 
a single hole A^ state (A ^ 1 —>>A ^ 2 transition). There¬ 
fore at zero bias, the observed acceptor pairs are in an 
A^~ charge state, and a bias V = aT 1 x 50 mV much 
less than —IV would be required to remove the last hole 
by tip-induced band bending. This is incompatible with 
U ~ 0.2 to 0.3 V of peaks observed in experiments. On 
the other hand, ionization energies of the two-hole states, 
A 9 —^A^ of coupled acceptors are calculated (Section [~C| 
to be very energetically similar to Ao ionization energies. 
Consequently, transitions populating ground and excited 
two-hole will occur for U >0, explaining how they (much 


d /ob 

Fit # 1 


Fit #2 



7oo 

S5T 

7oo 

S5T 

2.2 

free 

0.0159 

fixed, 0 

0.0293 

2.7 

free 

0.0126 

fixed, 0 

0.0509 

3.5 

free 

0.0391 

fixed, 0 

0.0638 


TABLE I. Comparison of sum of square errors when correla¬ 
tion parameter 7 00 is a free parameter in fitting model (Fit 
# 1), versus when 7 00 is fixed to zero (Fit #2). 


d(ae) 

2.2 

2.7 

3.5 

t (GHz) 

2900 

1700 

850 

Ms (GHz) 

< 2.5 

< 30 

< 10 

t/hTj: 

> 1200 

> 60 

> 85 


TABLE II. Estimated tunnel coupling 2 1 and environmen¬ 
tal coupling rates for each d. Conservative estimates are 
given for the molecular coupling, based on calculations for the 
[100] direction. Nevertheless, even the conservative molecu¬ 
lar coupling estimate for the [100] direction well exceed upper 
bounds on the environmental coupling, estimated from the 
dl/dU lineshape. 

like single-acceptor ionizing resonanceJ 29 ^) are observed in 
our experiments. 

2. Transport coupling to reservoir 

We have used the dl/dU lineshape to obtain an upper 
bound on the tunnel coupling T^ = T[ n + r ou t to the en¬ 
vironmental reservoirs - the largest lifetime broadening 
energy scalJ^ hT^ that can be included in the lineshape 
analysis of the data 29 . As necessary to observe molecu¬ 
lar states, we find for data in Fig. [3]4, [3^>, and[3p, 
to be considerably smaller than the inter-acceptor tun¬ 
nel coupling t estimated from the hybridization energies 
(Table [TT[). Tip height-dependent measurements (Section 
m confirmed the exponential dependence of the tunnel 
current on tip height, indicating the rate to the tip is the 
slower (dominant) rate. The dwell-time (and coupling 
r E ), determined by the fast rate, is therefore controlled 
by the coupling to the hole reservoir in the substrate, Tin 
in Fig. [6^1. 

B. Acceptor depth estimate 

We have estimated the depth below the silicon surface 
of the acceptor dopant pairs, presented in Fig. HM1P and 
[3p of the main text, to be approximately 0.5 nm, 0.9 nm, 
and 0.9 nm, respectively. The depth of isolated dopants 
in scanning tunneling spectroscopy is typically estimated 
fitting the spatial variation of a spectral feature such as 
a band edge, in a bias condition where the dopant is 
ionized 29 31 32 IH( to a dielectric screened Coulomb poten¬ 
tial. However, as illustrated in Fig. [6]E, tunneling from 
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FIG. 8. Dependence of the full-width half-maximum W of the 
spatial tunneling probability distribution on depth of acceptor 
beneath the surface, calculated by numerical diagonalization 
of the 6x6 spin-orbit coupled Kohn-Luttinger Hamiltonian. 


the reservoir to the acceptor pair is a transition from 
a one-hole state to a two-hole state. Consequently, for 
voltages U below the lowest voltage peak at U = U\ 
in Fig. and Fig. [6j3, the two-acceptor system is in 
Coulomb blockade with single-hole occupation, rather 
than being doubly ionized, such that fitting to dielectric- 
screened single-ion potentials is inappropriat e! 29 * 31 * 32 * 54 *. 
For U above U\ the two-acceptor system has a high 
probability of occupation by two holes since r ou t Tin. 
Hence, the method of fitting the band profile to dielec¬ 
tric screened ion potentials used in the above references 
is problematic. 

To estimate the depth of our acceptors, we use the 
well-known result that the spatial extent of the bound 
state measured by STS increases with increasing dopant 
depth (see e.g ., reference l55l). We compare the full-width 
at half-maximum W of the spatial tunneling distribu¬ 
tion fit in Section [D] to the same quantity calculated for 
subsurface acceptor-bound holes in the 6x6 spin-orbit 
coupled valence band. We provide this estimate for com¬ 
pleteness, since the fitting procedure described in Section 
[D]does not rely on prior knowledge of the depth or Bohr 
radius of the dopant. 

The full-width at half maximum W for the Is orbital 
density (j) s { r) ^ exp( —|r — ro|/a*) at z = 0 can be di¬ 
rectly solved giving W = (^a* log(|) — Zq) 2 — (zq) 2 . 

Evaluating this quantity we obtain W * 1.48 ± 0.12 nm, 
W = 2.25 =b 0.39 nm, and W = 2.19 =b 0.25 nm for the 
measured molecules with d/ap = 2.2, 2.7, and 3.5, re¬ 
spectively. Comparing these results to the full-width at 
half maximum W predicted for the orbitals using the nu¬ 
merical solution to Kohn-Luttinger 6x6 Hamiltonian 
presented in Fig. [8j the W extracted from experiments 
translate to depth estimates of ~ 0.53 nm, ~ 0.9 nm, and 
~ 0.9 nm, respectively. 


1. Correlation with lever arms 

When the acceptors are located deeper in the sample, 
the lever arms are anticipated to be smallei^ 9 . The rela¬ 
tive depths of the planes in which the acceptor pairs are 
estimated to be located agree with the extracted lever 
arms. We obtain a lever arm a = 0.0144 =b 0.0006 for 
the shallower d/ap = 2.2 pair. The other two d/ap = 2.7 
and d/a-Q = 3.5 pairs which are slightly deeper have lever 
arms a = 0.0088 d= 0.0010 and a = 0.0055 d= 0.0003, re¬ 
spectively. We note that fluctuations in the depth of the 
reservoir can also influence the lever arm. We anticipate 
that varying the dopant depths has a small effect on cor¬ 
relations in our experiments where the overall system is 
neutral. We have observed that the neutral level of near¬ 
surface acceptors (N = 1) remains bulk-lik^ 2 ® for depths 
0.5-2.0 nm, as expected, due to a competition between 
the confinement and the dielectric mismatch. Fixed ion¬ 
ization energy essentially fixes the Bohr radius, fixing t 
and U/t. We note however that the apparent spatial ex¬ 
tent of the wavefunction increases with increasing depth. 
This is not because of a change in the physical Bohr ra¬ 
dius. Rather, it occurs because STM probes the wave- 
function near the surface. 


C. Theory of interacting acceptors 

We employed a microscopic configuration interaction 
framework to calculate eigenstates of interacting holes 
bound to proximate acceptors, in a Kohn-Luttinger 6 x 
6 spin-orbit coupled representation (heavy holes, light 
holes, and split-off holes), and in a simplified single-band 
representation for reference. In both representations, the 
two-hole ground state at relatively small d/a b ^ 2 was 
found to be predominantly composed of the |e + e—) sin¬ 
glet of two even orbitals, where + and — denote +“3/2” 
and —“3/2” respectively for the valence band holes, and 
denote t and ^ respectively in the single-band approxima¬ 
tion. The probability amplitude of the |o + o—) singlet 
of odd orbitals was found to increase with increasing d, 
signaling interaction-driven generation of quantum corre¬ 
lations. The corresponding spatial tunneling probability 
distribution was predicted using the quasi-particle wave- 
function (QPWF) density appropriate for few-particle 
states with quantum correlation J 19 * 20 *, and are presented 
herein for the non-trivial 6x6 spin-orbit coupled case. 

1. Overview 

In the absence of spin-orbit coupling, wavefunctions 
of interacting spins can be classified into singlets and 
triplets 5 ®. For holes in the valence band, the atomic or¬ 
bitals have p-like symmetry, and the angular momentum 
L = 1 couples to the spin (intrinsic) angular momentum 
S = 1/2, such that the resulting Bloch states have defi¬ 
nite J and J z . While intrinsic angular momentum is not 
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a good quantum number, spin can nevertheless be gener¬ 
alized to hole pseudospin, z.e., Kramers doublets of time 
reversal symmetric linear combinations of Bloch stated 7 . 

Coupling of hole spins is well understood in the 
Kohn-Luttinger framework 5 ^®. It has been theoret¬ 
ically shown that coupled hole pseudospin interact to 
produce singlets and triplet^®, as for coupled electron 
spins. In the Kohn-Luttinger approach generalized va¬ 
lence band holes spins as written as spinors <^(r) = 
Ej, mj F j,mj( r )\ J ’ m j), where are the envelope 

functions corresponding to the Bloch states | J, raj). Due 
to spin-orbit coupling such spinors are ^ 90 % polarized 
(rather than 100 % polarized) into single Bloch compo¬ 
nents. This produces a small band-mixing effect when 
the pseudospins are coupled togetheii®®, that we take 
into account exactly in our numerics. 

2. Full-Configuration Interaction 

We employ the Kohn-Luttinger FCI approach devel¬ 
oped for valence band holes in reference [59] . Two- 


T E'ox'jdisx T Jvx K v x T ^ ^ d a 

ol/3 


hole eigenstates |T) are expressed as superpositions of 
antisymmetrized two-hole Slater determinants | a/3) = 

441°> = 2_1/2 (l a )i l/?> 2 - \P)i \ a )i) with probability 
amplitudes d a p, 


|M/) = ^d Q/3 |a/3). (1) 

a, (3 


In | a/3) above, cj creates a single-hole state | j) = cj |0), 
where |0) is the vacuum state. As discussed in detail be¬ 
low, the single-particle states | j) were chosen as the eigen¬ 
states of the non-interacting Hamiltonian FLo of a pair of 
acceptor ion potentials and a hard-wall interface with 
dielectric mismatch. The Hamiltonian for the two-hole 
system was taken as Ft = FLo (iq ) + FL o(r 2 )+ T4-h (r 1 , r 2 ) 
where Vh-h(i*i, r 2 ) is the Coulomb repulsion term for the 
holes. Substituting |T) above into FL\^) = F?|T) and 
taking the inner product with (z/A|, we obtain 


^a/3,ux) (Y vX,a.(3 Y u x,/3a') — -^Wa, (2) 


v\,ap = J dr l( lr 2 ^ J2 F; jiTO .(r 1 )F a) j im .(r 1 ))y h _ h (r 1 ,r 2 )( ^ (r 2 )), 

J,rrij J' ,m'. 


(3) 


for each \i/X), where Y ux ,a(3 = (^|i(A| 2 T4-h|^)i|^)2 is an 
off-diagonal (interaction) term evaluated in the spinor 
representation, J v ,x = Y u x,uX is the direct Coulomb in¬ 
teraction, and K v x = Y v x,Xv is a direct exchange interac¬ 
tion. Finite Y u x, a p terms can introduce correlations 59 by 
creating eigenstates that are admixtures of two-particle 
Slater determinants | a/3). 

The single particle Hamiltonian employed was FLo = 
^(k) T ldon,A T Wm,B T ^dmage,A T Mmage,B T Fw a ll taking 
into account ion potentials Wn,A/B of both acceptors, im¬ 
age charges Ki mage? A/B due to dielectric mismatclJ®, and 
an infinite (hard-wall) Wail potential a distance lap ^ 1 
nm away from the ion representing the semiconduc¬ 
tor/vacuum interface. For FL( k) we considered bot h the 
6x6 spin-orbit coupled Kohn-Luttinger valence band 61 62 
FLkl (k) , as well as a single-band, isotropic, parabolic dis¬ 
persion relation FL P ( k) = h 2 k 2 /2m* where m* is an effec¬ 
tive mass chosen to reproduce boron’s binding energ}E31. 
The latter allows us to calculate states of a simple “scaled 
hydrogen molecule” without complex spin-orbit coupling 
and band structure effects. The ion potentials were taken 
as dielectric-screened Coulomb potentials. For the 6x6 
scheme we included additional charge Sq = — O.le at the 
ion site as a crude core-correction to obtain the correct 


spatial extent ( r ) and ionization energy of the stated®. 
The interaction term Vh-h(ris, r 2 ) is a dielectric screened 
Coulomb interaction including an image charge distribu¬ 
tion produced by dielectric mismatch 6 ®. 

For the single band theory, the four lowest single¬ 
electron eigenstates of FLo were used to construct 6 
possible two-particle configurations, all of which were 
included for a full configuration interaction (FCI) ap¬ 
proach. We considered two even, spin-degenerate or¬ 
bitals |ecr) = c| a |0), and two odd, spin-degenerate or¬ 
bitals 10(7) = |0). 

For the 6x6 Kohn-Luttinger calculations, we use a ba¬ 
sis that reflects the four low-energy states of single accep¬ 
tors with s-like envelopes - the predominantly \mj\ = 3/2 
(“3/2]D state, and the predominantly \mj\ = 1/2 (“1/2”) 
state®. Near an interface, the “1/2” state experiences 
more confinement than the “3 /2” state, making the “3/2” 
state the ground statJ^ 0 67 EH > p or lone subsurface accep¬ 
tors ~ 1 nm from an Si[001]:H interface, we have typically 
measured ^1 — 2 meV splittings (see e.g ., lone accep¬ 
tor data in Fig. [7]), in good agreement with few meV 
values predicted by 6 x 6 Kohn-Luttinger. See reference 
( 30] for further details. In the two-acceptor potential FLo 
these hybridize into eight single-hole eigenstates, which 
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were used to construct 28 possible two-hole configura¬ 
tions, all of which were included in the state vector for 
a FCI approach. For the acceptor distances in the main 
text, four were predominantly composed of heavy holes 
(“3/2” states) and four were predominantly composed of 
light holes (“1/2” states). In both groups of four states, 
two are found within a Kramers doublet with even par¬ 
ity symmetry for the majority spin, and two are found 
within a Kramers doublet having odd parity symmetry 
for the majority spin! 58 59 . 

A numerical representation of the FCI matrix (Equa¬ 
tion [2| was directly diagonalized to obtain state vectors 
d and energies E. The FCI matrix was obtained by eval¬ 
uating the Coulomb matrix elements J u a, exchange ma¬ 
trix elements K v \, and the interaction terms Y u \^ by 
Monte-Carlo integration, using the VEGAS method for 
adaptive sampling and 10 7 iterations per integral. The 
basis of single particle eigenstates required for the inte¬ 
grals was obtained numerically using the finite difference 
scheme, for both the single-band and 6x6 representa¬ 
tions, on a real-space grid of 110 x 110 x 110 points with 
a discretization step size of 0.21 nm, using Luttinger pa¬ 
rameters for silicon’s valence bancP^. To normalize the 
measured and theoretical (6x6 and single-band) inter¬ 
acceptor distances, we choose the value ae = 1.3 nm 
reproducing the 44.4 meV binding energy of a Boron ac¬ 
ceptor in siliconP^ in a single-band approximation. 


3. Correlations in ground state 

In this section, the ground state of the FCI Hamil¬ 
tonian in Section |~C2] is discussed. Our FCI predicts a 
singlet ground state for both the spin-orbit coupled 6x6 
representation and for the isotropic single-band case. In a 
single-band representation (ignoring spin-orbit coupling), 
the ground state 

l*s> = (7ee4,t C e,i “ 7oo4 it <±)|0} (4) 

was obtained. With increasing d, dynamical interactions 
(T’s in Equation [2| resulted in increasing 7 00 and in¬ 
creasing quantum correlations C = 2 | 7 00 | 2 , as predicted 

elsewhere™!] 

T he gr ound state of interacting, spin-orbit coupled 
hole d 58 * 59 -! in the 6 x 6 Kohn-Luttinger representation con¬ 
tains contributions for which majority pseudospin have 
even parity, 

^ 1 ^) 5 ( 5 ) 

and contributions for which for majority pseudospin have 
odd parity, 

^ ^(o,m J ),(o,m' J ) C o,mj C o,m / J 1^) ’ (^) 

For small d/a b ^ 2, the dominant configuration (ap¬ 
prox. 90 %) in the ground state singlet was found to 


be d(e, 3 / 2 )(e,- 3 / 2 )j while for increasing d, the proba¬ 
bility amplitude d( 0? 3 / 2 )(o,- 3 / 2 ) was found to increase. 
From this, two conclusions can be drawn: (1) the in¬ 
crease in the probability amplitude d( 0 j 3 / 2 )(o,- 3 / 2 ) with 
increasing d signals the emergence of hole-hole quantum 
correlation^, and (2) the predominance of “3/2” pseu- 
dospinspin orbitals signals that the band-mixing effect 
is weak for the values of d in experiments. Contribu- 
tions and d(o,mj),(o.m'j) with "1 j = ±3/2 

and m'j = ±3/2 are much larger than the other “band¬ 
mixing” terms with mj = ±“l/2’ or m'j = ±“1/2”. 


4- Quasi-particle wavefunction density 


In this section, we discuss theoretical predictions of 
the correlation coefficient C in Fig. [4}A, as well as theo¬ 
retical calculation of the quantity measured for spatially 
resolved single-hole tunneling into mul ti-pa rticle states 
with correlations, the QPWF densit y* 20 * 22 ! The results 
of Section C 3 are considered, for both the single-band 
and spin-orbit coupled 6x6 Kohn-Luttinger representa¬ 
tions. 

As described in Section 0 the tunneling probability 
density measured in experiments is determined by tun¬ 
neling of a hole from the acceptor molecule to the tip, 
leaving behind a single-hole eigenstate | i) on the acceptor 
molecule. The corresponding QPWF density Po Ut (r) oc 
\Mis(r)\ 2 is obtained^, where M^(r) = (i|^(r)|^), 
*(r) = E; c j4>j{ r) is the field operator, r = (x,y,zo) 
is position of the tip apex orbital, and \^s) is the 
ground state. For the single-hole transport processes 
probed, it is necessary to sum over the final (single¬ 
hole) states | i) of the molecule after hole tunneling to 
the tip. The total tunneling probability density is then 

r 0 ut(r) = Ei Uut(r)- 


For the single band case, evaluating the 
QPWF using the field operator gives r out (r) oc 
[l7ee| 2 |<Mr)| 2 + |7oo| 2 |0o(r)| 2 ]- Results for C = 2|7 GO | 2 
are plotted in the main text (Fig. Eh dashed line). 

For the 6x6 spin-orbit coupled representation, the 
quasi-particle wavefunction densities are also readily 
obtained, such that r out (r) = |M^( r )| 2 , where 

M iS (r) = E a/ 3 d a ^ a (r)S ij3 - </>b{r)8 ia ]. Evaluating this 
expression results in a QPWF density that is a sum of 
even and odd QPWF contributions for all d and both 
orientations (100) and (110), as expected considering the 
generalized parity of holes 58 ! Results for C, again defined 
as twice the probability density of odd orbitals, are plot¬ 
ted in the main text (Fig. coloured lines) for both 
(100) and (110) orientations. 

The predicted spatial tunneling probability density is 
shown 1 nm above the ion in Fig. [9] for an acceptor 
molecule with d = 3.8aB* The predicted correlation co¬ 
efficient is C = 0.6, similar to C = 0.78±0.16 extracted 
from measurements of the d = 3.5aB acceptor molecule, 
in Fig. [4jA of the main text. Contributions to the spatial 
maps (line profiles) from even and odd QPWF densities 










12 



x(nm) 


x(nm) 


x(nm) 


FIG. 9. A. Predicted spatial tunneling probability distribution for even quasiparticle wavefunctions using 6x6 Kohn-Luttinger. 
B. Same as A, but for odd quasiparticle wavefunctions. C. Total tunneling probability distribution, a sum of A and B. D. 
Profile of spatial tunneling probability distribution for even quasiparticle wavefunctions using 6x6 Kohn-Luttinger. E. Same 
as D, but for odd quasiparticles wavefunctions. F. Total tunneling probability distribution, a sum of of D and E. Scale bar: 1 
nm 


are separately plotted in Fig. [9]A (Fig. [9jD) and Fig. [9p 
(Fig. |9p), respectively. Their sum, the total tunneling 
probability density (line profile), is plotted in Fig. [9p 
(Fig. [9p). Upper and lower grey lines in Fig. [9p denote 
normalized QPWF densities for the even QPWF and odd 
QPWF components, respectively, which in this case di¬ 
rectly correspond to the profiles in Fig. [9p and [9p, re¬ 
spectively. 


D. Fitting of measured spatial tunnelling 
probabilities 

The experimentally obtained spatial tunneling prob¬ 
abilities for the two-hole ground state of acceptors, in 
Fig- [3]A., [3j3, and[3p of the main text, were fit assuming 
even and odd linear combinations of atomic orbitals with 
s-like envelopes for 0 e (r) and 0 o (r). This data was ob¬ 
tained using a tip height established by constant current 
imaging conditions, at a bias where the topography was 
nominally flat apart from dimer row corrugations. 

In this section we describe the fits and show that the 
probability density of an s-like wavefunction describes the 
measured probability density of the ground state and first 
excited state of an isolated acceptor, in agreement with 
what is expected for the s-like “3/2” pseudospin ground 
state and s-like “1/2” pseudospin excited state. 

Spatially resolved tunneling spectra were measured 
for isolated acceptors, consistently revealing two dl/dU 
peaks associated with the s-like J — 3/2 ground state 
manifold. The ^ 2 meV splitting of the isolated accep¬ 
tor in Fig. [7] was found to be in good agreement with 
single-particle eigenstates calculated from our numerical 
Kohn Luttinger solver for depths determined as in previ¬ 
ous work^l, that is, by a fitting the profile to t he ionized 
acceptor’s perturbation to the band edge 32 54 l The spa¬ 


tial tunneling probability distributions for the lowest en¬ 
ergy state and first excited state of an isolated acceptor 
were obtained by integrating the two low energy d //dU 
peaks. Results are shown in Fig. [T()] A and Fig. [Top for 
the single acceptor in Fig. [7| 

The dimer lattice structure of the hydrogen-terminated 
2x1 surface is evident in the states Fig.[To|4 and Fig.flOp, 
running along the [110] direction indicated. The dimer 
and lattice frequencies were Fourier filtered out of the 
imagd^, after which only the characteristic envelope of 
the probability density of the states remain. Results for 
the ground state and first excited state are presented in 
Fig. [Top and pTp, respectively. The tunneling probabil¬ 
ity density reflects the hole density, and the appearance of 
two s-like envelopes is in agreement with expectations of 
the s-like envelopes for the “3/2” and “1/2” state^^. 
We observe that the measured low-frequency envelope 
for the probability density is closer to isotropic for the 
ground state, and spatially elongated along the dimer 
direction for the excited state. 

The probability density was fit along a line in the 
plane (z = 0) of the surface of the form |0 s (r)| 2 , where 
4> s (r) ~ exp(— |r — r 0 |/a*) is a Is orbital envelope func¬ 
tion. In this expression, r = (x, y , z) is the tip position, 
r 0 = (xo, 2/cb ^o*)> x o and 2/0 are centre coordinates of 
the orbital, Zq is an effective depth, and a* is an effec¬ 
tive Bohr radius. The profile of the measured ground 
state (solid squares) and least-square fit (solid line) are 
shown for a profile along the direction parallel (black) 
and perpendicular (green) to the dimer row in Fig. [Top, 
demonstrating excellent agreement. The same quanti¬ 
ties are plotted for the first excited state in Fig. [lOp. A 
more noticeable anisotropy of the excited “1/2” state was 
found. 

Motivated by the theoretical observation that the 
QPWF density partitions into contributions from even 
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FIG. 10. A. Integrated dl/dU for single acceptor ground 
state. B. Integrated dl/dU for first excited state, 1.6 d= 0.4 
meV above the ground state. C. Same as A with lattice fre¬ 
quencies filtered out. D. Same as B with lattice frequencies 
filtered out. E,F. Line profile of relatively tunneling probabil¬ 
ity parallel (black squares) and perpendicular (green squares, 
offset for clarity) to the dimer, and fit of profile parallel (black 
line) and perpendicular (green line, offset for clarity) to dimer 
direction. F. Same as E, but for first excited state. Scale bar: 
1 nm 


and odd states, we categorize the measured QPWF den¬ 
sity into contributions from even and odd QPWFs, that 
is, linear combinations of atomic orbitals with even and 
odd parity, respectively. For a distance R between accep¬ 
tors, we use 0 e (r) = (l// e (i?))[0 a (r-R/2)+0 a (r+R/2)] 
and cf) 0 ( r) = (l/J o (R))[0 s (r - R/2) - 0 s (r + R/2)], 
where J e (r) = 1 + (1 + r + r 2 /3)exp(—r) and J 0 (r) = 
1 — (1 + r + r 2 /3)exp(—r) in atomic units. In the fits 
presented in the main text, a*, Zq, R, and |7 ee | 2 are free 
parameters, while 1 = |7 00 | 2 + l7ee| 2 is simply a con¬ 
sequence of normalization. Effective Bohr radii a* ^ 1 
nm and depths ~ 1 nm were obtained from the least 
squares fits. In our fits the values for a* and Zq define a 
full-width at half maximum W for the spatial tunneling 
probability of each acceptor. In Section m the actual 
depth zq of the acceptors was estimated using W derived 
from experimentally fit values a*, Zq. 
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FIG. 11. A. Spatial measurement of ground state for d = 
2.2as- B. Same as A, but for first excited state. C. Spatial 
measurement of ground state for d = 2.7an. D. Same as C, 
but for first excited state. (110) crystal directions directions 
are as indicated. Scale bar: 1 nm 


1 . Error analysis in fit 
The sum of square errors 

= 57 r (*i) - r fl t(*i))7(r(*i)) 2 (7) 

Xi 

was used to determine the importance of the Coulomb 
correlation parameter 7 00 in the fits. The first column 
(Fig ffl) of Table |l| gives the results for the SSe when 
|7oo is a free parameter, while the second column (Fig 
# 2) of Table [i] gives the results for SSe when | 7 00 is 
artificially forced to zero. Here we see the fit is consid¬ 
erably improved by including the Coulomb correlation 
parameter in the fit. 


2. Independence of fit results on tip height 

Due to the very small lever arm a ~ 10 —2 , the data 
is acquired very near flat-bancP^ and as such we do 
not observe the so-called “ionization parabolas” (ob¬ 
served elsew here and associated with tip-induced band 
bending 2 33 ^) for single B:Si acceptors 3 ^, or for B:Si ac¬ 
ceptor pairs (Fig. lb, main text). 

To provide an independent check that the tip does 
not influence the results of the QPWF measurement, we 
present measurements on the d/a b = 2.2 and d/a-g = 3.5 
pairs from the main text at a tip height +60 pm above the 
position where the data was taken in the main text. The 
normalized tunneling probability is shown in Fig. [l2] A 
for d/a-B = 2.2 and in Fig. ["L2^ 3 for d/a b = 3.5. Line 
profiles of the normalized tunneling probability (squares, 




















14 



x(nm) 


x(nm) 


FIG. 12. A. Experimentally measured, normalized tunneling probability T to tip, for d — 2.2as GS. Arrows denote 110 crystal 
directions. B. Same as (A), for d = 3.5aB. C. Normalized experimental line profile (coloured squares) of r(x) for d m 2.2 clb 
and least-squares fit (coloured line) to correlated singlet model. Upper and lower grey lines are line profiles of densities |0 e (r)| 2 
and 1(r) | 2 respectively, obtained from least square fits. D. Same as (C), for d = 3.5 as- Scale bar: 1 nm 


Fig. [~L2] G and Fig. [~L2p), are in excellent agreement with 
the same QPWF model (solid red line, Fig. [TTp and 
Fig. |l2p) used in the main text. Reference curves are 
shown for the totally uncorrelated state ( 7 00 = 0) and 
maximally correlated state (7 00 /7ee = 1). 

The extracted results for the correlations are 17 00 1 2 = 
0.15 ± 0.06 for d/a B = 2.2 and | 7 00 | 2 = 0.40 ± 0.09 for 
d/a B = 3.5. These results for +60 pm tip heights are 
within the experimental errors of the results for the data 
in the main text, which are 17 00 1 2 = 0.12+0.06 for d/a B = 
2.2 and |7 00 | 2 = 0.39 ± 0.08 for d/a B = 3.5. 


3. Large distance limit of fitting model 

For any given experimental signal-to-noise ratio, there 
always exists a large enough d/a B where it is impos¬ 
sible to distinguish the existence of correlations in the 
ground state QPWF. This follows from the fact that 
our model determines correlations based on the con¬ 
trast between orbitals |0 e (r)| 2 and |0 o (r)| 2 . However, 
since |<^ e (r)| 2 - |</>o(r)| 2 ~ 2</> 7 i(r)^ B (r), where 4> A /b( r) 
are atomic orbitals at sites A and F>, and the product 
20^(r)0s(r) is exponentially suppressed with increasing 
d/a B . The results of Table IT] show that the Coulomb cor¬ 
relations embodied by | 7 o 0 Tcan be reliably detected, for 
the signal-to-noise ratio of our experiment. 


4- Spatial measurements of excited states 

We present the ground and excited states of the accep¬ 
tor pairs for d/a B = 2.2 and 2.7 in Fig. [TT] Compared 
with the ground states for d/a B = 2.2 and d/a B = 2.7 in 


Fig. [TTT A and Fig. [UpS respectively, the excited states 
for d = 2.2 a B (Fig. W) and d = 2.7 a B (Fig. ET) 
are slightly spatially elongated along the 110 direction. 
This is consistent with the expected QPWF density in 
the interpretation of Fig. [5] (main text), that the excited 
state orbitals have one hole occupying a “1/2” pseudospin 
state and the other occupying a “3/2” pseudospin state. 
These states display slightly different spatial anisotropies 
(Fig. [Top and Fig. [Top). ' 

While we have measured the spatial tunneling prob¬ 
abilities for the excited states, measuring excited-state 
quasi-particle wavefunctions requires slow tunnel-in from 
the tip into the dopant system (not slow tunnel-out to 
the tip). Slow tunnel-in from the tip is required to make 
the current limiting slow tunnel rates additive for excited 
states (rather than competitive with the ground state as 
in the case herein), in the presence of Coulomb blockade. 
Some details are given in reference [40j 


E. Tip-height dependent spectra 


In this section we discuss our confirmation of our 
model for the thermally broadened sequential tunneling 
in Fig.[6p, by detailed investigation of the tip-height de¬ 
pendence of the curr ent, in agreement with our results on 
Arsenic donors in Si 3140 . This is shown in Fig. 13 A for 
tip heights that sequentially increase by 15 pm per curve, 
above a single B:Si acceptor presenting the unique topo¬ 
graphical signature discussed elsewhere 29 30 . The data 
shows broadened current steps centred at U « 0.25 V 
and U w 0.5 V that have an exponential dependence 
on the tip height, and which fit very well the lineshape 
of purely thermally broadened current resonance of a 
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FIG. 13. A. Tip height dependence of current I vs sample bias 
U for single acceptor. Two distinct thermally broadened steps 
corresponding to ±3/2 states and ±1/2 states are observed 30 . 
B. Peak current Io reduces exponentially as a function of tip 
height. C. Energy splitting to the ±1/2 excited state 


F. Hubbard model vs. Molecular orbitals 

The results for the two-site Hubbard model and its 
relationship to the molecular orbital model in Fig. [2^3 
of the main text are straightforward to derive. The co¬ 
efficients 7 c and 71 of the singlet ground state \^s) = 
7c(|m) - II; t» +7i(ltl;> + |;TI» (Fig.|t, main text) 
are readily found by direct diagonalization to be q c = 
a/V 2 ± a 2 and 7 * = l/\/2~± ~2a? where a = {—U/At ± 

yi + WJW)~ 1/2 - 

The dependence of 7 ee and 7 00 on U/t in Fig. [2p can 
be found by a simple transformation. The orthonormal 
localized states of the Hubbard dimer model created by 


operators c^ Aa and c^ Ba are 


C Aa = W ~ 1/2 ( C L 
C Ba = W ~ 1/2 ( C L 

- gcL) and 

- g<i*) 


respect ivel}!^. Here cI aG and create (non-orthogonal) 
atomic orbitals with spin < 7 , w = 1 — 2 Sg + g 2 , g = 
(1 — \/\ — S 2 )/S , and S = (a\b) is the overlap between 
the atomic orbitals. 

Re-writing the even and odd eigenstates of the single¬ 
particle potential and equation the two singlet ground 
states we find 7 ee = 7 C ± 7 i and 7 00 = 7 c — 7 i • These rela¬ 
tionships are used to map the Hubbard model solution to 
the molecular orbital model in Fig. [2^3 of the main text. 
This mapping is independent of normalization constants 
re, g , and S. 


dopant ! 29 31 4 Q l 53 i The peaks shown here correspond to 
the ±3/2 ground state and ±1/2 first excited state dis¬ 
cussed in Section m and reference [30l 


Qualitatively, the exponential dependence on tip 
height verifies the assumption that T[ n r out , as ex¬ 
pected for the large vacuum tunneling barrier, and as 
required to obtain the results in Fig. [3] of the main text. 
We performed least-squares fitting to extract the height 
Io of the thermally broadened current steps. These Iq 
are plotted as a function of tip height z - Z 0 in Fig. [13 }d 
for the ±3/2 ground state (red curve) and the ±1/2 
excited state (blue curve). The data are in excellent 
agreement with a z-dependence exp (—2 k(z — zq)) where 
k = (1.15 ± 0.05) x IO 10 m -1 . This value is in good 
agreement with the expected value for k = \/ 2 mo$B/S 
for the ~ 5 eV barrier expected for shallow valence band 
acceptor in silicon. Moreover, we find that the lever arm 
(not shown) and energy splitting (Fig. [T3]G) between the 
ground and excited state does is essentially independent 
of tip height. This provides further verification of the 
sequential tunneling model, and that the tip does not 
strongly influence the acceptor-bound states. 


G. Degree of entanglement 


For fermions, entanglement is defined in terms of Slater 
decompositions 23 rather than the Schmidt decomposi¬ 
tions of distinguishable particles, and can be expressed 
independent of basis by a degree of entanglement. Here 
we follow derivation in reference [36] for entanglement 
to describe correlations in a two-site Fermi-Hubbard sys¬ 
tem. For two particles that obey Fermi statistics, 


w = X±±> i i & > 2 


( 10 ) 


a,b 


is characterized by an anti-symmetric matrix c o a ,b, that 
is, a Slater decomposition. Transformed into a block di¬ 
agonalized form diag[Zo, Zi ,..., Zn\ through a unitary ro¬ 
tation of the single particle states, where 


Zi = 




(ii) 


the number of nonzero Zi is called the Slater rank, and 
if the Slater rank is 1 , the quantum correlation of the 
state is zercP^. It has been shown 36 that the z\ are the 
eigenvalues of the basis independent quantity cohu;, and 
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based on this observation, the degree of entanglement S 
of the two particles was defined as 


S = - J2 Z i l °&2 Z i- ( 12 ) 


glets \S) = (7ee4 >t 4 4 - 7 o<4 it 4 )4 _)|0>, 


^ 0 7ee 0 0 ^ 

7 ee 0 0 0 

0 0 0 7oo 

\ 0 0 7oo o / 


(13) 


It directly follows that the non-zero eigenvalues zf of lj^lu 
are 17 ee | 2 and |7oo| 2 - Then, the degree of entanglement 
for the state |£) is 


For a superposition of “even/even” and “odd/odd” sin- 


s = -|7ee| 2 log 2 |7ee| 2 “ |7oo| 2 log 2 |7oo| 2 - (14) 
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